version 19
drop _all

* change the root folder according to your computer directory
global root = "D:\WorkingPaper-Series\Township_Light\Replication_Files"  

global dofiles = "$root\dofile"     
global working_data = "$root\workdata"
global tables = "$root\table"
global figures = "$root\figure"

cd "$working_data"


******************
*****Pre Set******
******************

global geo_ctrl = "topography_t* slope_t* roughness_t*"
global wth_ctrl = "gaez_t* precipitation_t* sunlight_t*"
global dist_ctrl = "county_dist_t* dist_hz_t* dist_coast_t*"

global full_ctrl = "$geo_ctrl $wth_ctrl $dist_ctrl"

global geo_ctrl_1 = "topography_t_1 slope_t_1 roughness_t_1"
global wth_ctrl_1 = "gaez_t_1 precipitation_t_1 sunlight_t_1"
global dist_ctrl_1 = "county_dist_t_1 dist_hz_t_1 dist_coast_t_1"

global full_ctrl_1 = "$geo_ctrl_1 $wth_ctrl_1 $dist_ctrl_1"

global geo_ctrl_2 = "topography_t_2 slope_t_2 roughness_t_2"
global wth_ctrl_2 = "gaez_t_2 precipitation_t_2 sunlight_t_2"
global dist_ctrl_2 = "county_dist_t_2 dist_hz_t_2 dist_coast_t_2"

global full_ctrl_2 = "$geo_ctrl_2 $wth_ctrl_2 $dist_ctrl_2"

global matchvar = "county_dist topography slope gaez sunlight dist_hz dist_coast"


global pre_ctrl = "light_t road_t pop_t"


********************
***** Figure 3 *****
********************

* TWFE OLS

use "township_panel_final.dta", clear

gen g_5p = period < -5
    forvalues k = 5(-1)1 {
       gen g_`k' = (period == -`k')
    }
    
    forvalues k = 0/6 {
         gen g`k' = (period == `k')
    }

replace g_1 = 0

reghdfe log_light g_5p-g6 $full_ctrl treat_trend, a(county_year idcode) cl(countycode)

forvalues i = 1/12 {
    local m_`i' = e(b)[1,`i']
}

mat input mat1_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10', `m_11', `m_12')

mata st_matrix("V",sqrt(diagonal(st_matrix("e(V)"))))

mat V = V'

forvalues i = 1/12 {
    local m_`i' = V[1,`i']
}

mat input mat2_fe = (`m_1',`m_2', `m_3',`m_4',`m_5',`m_6',`m_7', `m_8', `m_9', `m_10', `m_11', `m_12')


mat twfe = [mat1_fe\mat2_fe]
mat twfe = twfe'
svmat twfe

keep if twfe1 != .
keep twfe1 twfe2
rename (twfe1 twfe2) (coef std_err)

export delimited using "twfe_est.csv", replace

* dCDH 2020

use "township_panel_final.dta", clear

cap drop light_resid
reghdfe log_light $full_ctrl treat_trend, a(county_year idcode) resid
predict light_resid, resid

did_multiplegt light_resid first_year year village, robust_dynamic breps(500) cluster(countycode) dynamic(6) placebo(5)

mat mat1_dcdh = e(estimates)'
mat mat2_dcdh = J(1,12,0)
forvalues i = 1 (1) 12 {
    mat mat2_dcdh[1,`i'] = sqrt(e(variances)[`i',1])
}

mat dcdh_est = [mat1_dcdh\mat2_dcdh]
mat dcdh_est = dcdh_est'
svmat dcdh_est

keep if dcdh_est1 != .
keep dcdh_est1 dcdh_est2
rename (dcdh_est1 dcdh_est2) (coef std_err)

export delimited using "dcdh_est.csv", replace

* SA 2020

use "township_panel_final.dta", clear

gen g_5p = period < -5
    forvalues k = 5(-1)1 {
       gen g_`k' = period == -`k'
    }
    
    forvalues k = 0/6 {
         gen g`k' = period == `k'
    }


drop g_1

gen never_treat = (treat==0)
eventstudyinteract log_light g_5p-g6,  cohort(first_year) covariates($full_ctrl treat_trend) control_cohort(never_treat) absorb(idcode county_year) vce(cl countycode)

mat mat1_sa = e(b_iw)[1,1..12]'
mata st_matrix("V",diagonal(st_matrix("e(V_iw)")))

mat V = V[1..12,1]

mat mat2_sa = J(1,12,0)

forvalues i = 1 (1) 12 {
    mat mat2_sa[1,`i'] = sqrt(V[`i',1])
}

mat mat2_sa = mat2_sa'

mat sa_est = [mat1_sa,mat2_sa]
svmat sa_est

keep if sa_est1 != .
keep sa_est1 sa_est2
rename (sa_est1 sa_est2) (coef std_err)

export delimited using "sa_est.csv", replace

* CS 2021

use "township_panel_final.dta", clear

reghdfe log_light $full_ctrl treat_trend, a(idcode county_year) resid
predict light_resid, resid

csdid light_resid , drimp i(idcode) t(year) g(first_year) agg(event) cl(countycode)

mat mat1_cs = e(b)[1,7..18]'
mata st_matrix("V",diagonal(st_matrix("e(V)")))
mat V = V[7..18,1]

mat mat2_cs = J(1,12,0)

forvalues i = 1 (1) 12 {
    mat mat2_cs[1,`i'] = sqrt(V[`i',1])
}

mat mat2_cs = mat2_cs'

mat cs_est = [mat1_cs,mat2_cs]
svmat cs_est

keep if cs_est1 != .
keep cs_est1 cs_est2
rename (cs_est1 cs_est2) (coef std_err)

export delimited using "cs_est.csv", replace
